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The BSSN (Baumgarte-Shapiro-Shibata-Nakamura) formulation of the Einstein evolution equa- 
tions is written in spherical symmetry. These equations can be used to address a number of technical 
and conceptual issues in numerical relativity in the context of a single Schwarzschild black hole. One 
of the benefits of spherical symmetry is that the numerical grid points can be tracked on a Kruskal- 
Szekeres diagram. Boundary conditions suitable for puncture evolution of a Schwarzschild black 
hole are presented. Several results are shown for puncture evolution using a fourth-order finite 
difference implementation of the equations. 



I. INTRODUCTION 

The Baumgarte-Shapiro-Shibata-Nakamura [TJ |2] 
formulation of the Einstein equations has been in 
widespread use in the numerical relativity community for 
a number of years. BSSN is the formulation used with the 
successful moving puncture technique for black hole sim- 
ulations [T, T . Currently several groups are using BSSN 
with moving punctures to model the gravitational wave 
signal from binary black hole mergers. (See, for example, 
Refs. laiilTllHliailllllllll].) 

There are still many open issues in numerical rela- 
tivity, both technical and conceptual. One would like 
to explore these issues in an efficient manner. Current 
three-dimensional codes can take days or longer on a 
multi-processor computer cluster to carry out each sim- 
ulation. By assuming spherical symmetry we can reduce 
the run time to minutes on a single processor. Thus there 
is an enormous practical advantage in assuming spherical 
symmetry. Of course, with a spherically symmetric (one- 
dimensional) code we are limited to issues that can be 
addressed in the context of a single Schwarzschild black 
hole. Furthermore there is no guarantee that the results 
one finds in spherical symmetry will transfer to the three- 
dimensional setting — in three dimensions there is a much 
wider range of possible problems that can arise. On the 
other hand, if a particular idea cannot be applied suc- 
cessfully in one dimension, it is unlikely to work in three. 

This article describes the one-dimensional code used in 
Refs. [inillllllS] to explore the evolution of Schwarzschild 
black holes. Such a code could also be used to examine 
various slicing and shift conditions, to explore different 
types of boundary conditions, and to experiment with 
constrained and unconstrained evolution schemes. Valu- 
able insights into these and other problems can be gained 
with the help of a one-dimensional code. 

Many questions that one can address with a one- 
dimensional code are independent of the formulation of 
the Einstein evolution equations. But it is reasonable 
to expect that the behavior of a one-dimensional code 
based on the BSSN formulation will be closer to existing 
three-dimensional BSSN codes. This is important when 
considering issues related to stability, constraint viola- 
tions, etc. Also note that the F-driver shift condition, 
currently used for puncture evolution studies of binary 



black hole coalescence, makes explicit use of the confor- 
mal connection functions F". The conformal connection 
functions are basic variables in the BSSN formulation. 

There is a technical obstacle to writing the BSSN for- 
mulation of Einstein's equations in spherical symmetry. 
More generally, there is a technical obstacle to writing 
BSSN in spherical coordinates. The obstacle is this: The 
BSSN equations assume that the conformal metric has 
determinant equal to one. This condition is not gener- 
ally covariant. In practice this is a problem because we 
typically use initial data that is conformally flat, and 
we naturally choose the initial conformal metric to be 
flat. But in spherical coordinates the flat metric does 
not have unit determinant. This obstacle can be over- 
come by generalizing the BSSN equations to allow for a 
conformal metric with non-unit determinant. The way 
to do so was presented in Ref. |16j . and is reviewed in 
Sec. 2. 

The generalized BSSN equations (GBSSN) are invari- 
ant under conformal transformations [16j . This invari- 
ance is broken when one chooses an evolution equation 
for the determinant of the conformal metric. There are 
two natural choices that yield the "Lagrangian case" and 
the "Eulerian case". In Sec. 3 I show how the GBSSN 
equations are reduced to spherical symmetry for each of 
these cases. The level of hyperbolicity for these equations 
is discussed in Sec. 4. 

In this work I focus on puncture evolution for a 
Schwarzschild black hole. The puncture is located at the 
coordinate origin. In spherical coordinates, the origin is 
also a coordinate singularity. Clearly the boundary con- 
ditions at the origin play an important role. In Sec. 5 
I describe sets of boundary conditions that work well in 
practice. In Sec. 6 I show how the numerical grid points 
can be mapped to a Kruskal diagram. This allows us 
to visualize the spacelike slice as it evolves, and track 
the distribution of grid points relative to the horizon and 
physical singularity. Section 7 presents some results for 
puncture evolution in one dimension. These results use 
the standard "l-|-log" slicing condition 

dta = fl^daa - 2aK , (1) 

or the related condition obtained by dropping the ad- 
vection term f3'^daa. Here, a is the lapse function, 
is the shift vector, and K is the trace of the extrinsic 
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curvature. For the shift vector, the examples use either 
vanishing shift or the F-driver condition 



(2a) 



dtB" = dtT" + p^dcB" - P'^dcT" - 7/B° . (2b) 

Variants of the F-driver condition are obtained by drop- 
ping one or more of the advection terms f3'^dcP"', P'^dcB"', 
or P^dcT" [7l[17|. 



II. GENERALIZED BSSN EQUATIONS 

The BSSN equations can be generahzed to allow for the 
possibility that the determinant of the conformal metric 
differs from unity |16j . The conformal metric gat and the 
trace-free part of the extrinsic curvature Aat are defined 

by 



gab = e'^'^gab , 



K, 



ah 



e'^ ( Aab + ^9abK 



(3a) 
(3b) 



where gab and Kab are the physical spatial metric and ex- 
trinsic curvature. The variable is the conformal factor 
and K = g°'''Kab is the trace of the physical extrinsic cur- 
vature. The conformal connection functions are defined 

by 
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^ g'^Tl = - — d, 
V9 
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(4) 



where F^^ are the Christoffel symbols built from the con- 
formal metric. The evolution equation for F° is obtained 
by differentiating this definition and using the momen- 
tum constraint. The BSSN variables are 0, gab, Aab, K, 
and F". 

In vacuum, the generalized BSSN evolution equations 
(the GBSSN equations) are^ 



d±cf> 

d±gab 

d±Aab 



-^9j_ln.g- ^aK , 
l^gabdj_lng - 2aAab , 
^AabdA_lng - 2aAacAl + aAabK 



TF 



-aK^ + aAabA'''' - D'^Daa , 



(5a) 
(5b) 

(5c) 
(5d) 
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2a 



,(5e) 



^ Equation jsj;) corrects a sign error in Ref. |16l . 



where Rab is the physical Ricci tensor. Da is the physi- 
cal covariant derivative, and Da is the conformal covari- 
ant derivative. The shift vector /?" is buried in the time 
derivative operator, d± = dt — Cp. Also note that TF 
stands for the trace-free part of the expression in brack- 
ets. 

The variables 0, gab, Aab, and K are defined as tensors 
with no density weights. The transformation rules for 
the conformal connection functions F° are determined 
by the transformation rules for the Christoffel symbols 
F^^ and the definition Eq. (|4|. The GBSSN equations 
do not assume g = 1, and we are free to choose how g 
evolves in time. Two natural choices are the "Eulerian 
condition" d±^ In g = and the "Lagrangian condition" 
(9(ln g)/9t = 0, which implies d±\ng = —2Da(3°'. 

In the traditional BSSN formulation, the restriction 
g — 1 appears as an extra constraint. In most numerical 
codes this constraint is enforced by replacing gab with 
9ab/9^^^ at the end of each time step. There is no clear 
equivalence between the traditional BSSN formulation 
and GBSSN with either the Eulerian (GBSSN-E) or La- 
grangian (GBSSN-L) conditions. In particular, observe 
that there are 16 field variables evolved by the GBSSN 
equations (namely 0, the six components of gab, the five 
components of Aab, K, and the three components of F") 
but only 15 variables in the traditional BSSN system. 
One consequence of this difference is a mismatch in the 
number of characteristic fields for GBSSN and traditional 
BSSN. In spite of this difference one might guess that GB- 
SSN with the Lagrangian condition dg/dt = is closely 
related to traditional BSSN. With GBSSN-L the deter- 
minant g, which can be set to 1 by initial conditions, 
remains constant throughout the evolution. Indeed, a 
qualitative comparison between the results presented in 
Section VII and results obtained in three dimensions with 
traditional BSSN shows that GBSSN-L and traditional 
BSSN yield very similar behavior for puncture evolution 
of a single Schwarzschild black hole. For example, com- 
pare Figs. 7-10 of this paper with Fig. 1 of Ref. [7]. 

Although GBSSN-L and traditional BSSN are similar 
in some respects, GBSSN-L contains an extra character- 
istic field. That field travels along the integral curves of 
the time flow vector field dt- Thus, this field has a char- 
acteristic speed that, when measured with respect to ob- 
servers at rest in the spacelike hypersurfaces, depends on 
the shift vector. This is not necessarily bad, but it is at 
least unusual and unphysical. For a typical formulation 
of the Einstein equations, the characteristic speeds are 
independent of the shift vector as long as any dynamical 
gauge conditions are expressed in terms of the normal 
derivative operator d± — dt ~ Cp. GBSSN-L breaks this 
pattern. On the other hand, for GBSSN-E, the extra 
characteristic is a zero-speed mode. As a result, the Eu- 
lerian condition yields a more simple and in some sense 
more physical version of the GBSSN system. 

Finally, note that the variable can be replaced with 
a new variable x = e"*"^. It is straightforward to rewrite 
the GBSSN equations in terms of x- Note that x, like 
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(f>, is a scalar with no density weight. The variable 
was used in Ref. [4], while the authors of Ref. [3j used 
X- Reference [7] includes comparisons between the two 
choices for black hole puncture evolution. 



and for the conformal connection functions, 



= I -cose/ {gee 9) 




(8) 



III. GBSSN REDUCED BY SPHERICAL 
SYMMETRY 

Let the spatial coordinates be denoted by r, 9, and Lp. 
The GBSSN equations are reduced to spherical symmetry 
with the following ansatz for the metric: 



grr 
gab= \ gee 

gee sin^ 9 



(6) 



The dynamical variables are now 4> or x, gm gee, A^r, K, 
and r*". They are all functions of the spatial coordinate 
r and time t. 

For the rest of this paper I will use the variable x rather 
than (f). I also use d±h\g = —2vDa/3"^ where v = gives 
the Eulerian condition and v — I gives the Lagrangian 
condition. In spherical symmetry, the GBSSN equations 
are as follows: 



The ansatz for the trace-free part of the extrinsic curva- 
ture is 



1 

Anh ^ \ -gee/{2grr) 

-gee&m^ 9/{2gr 



(7) 



= ^ - - - \vrx + , (9a) 

6 6grr ogee o 

a o /I ^ I , nr I '^grr'^P gee , n ar' ^ ar' /•nu\ 

Otgrr = -2A„a - -VP grr + P 5rr 5 1" ^grrP " l^grrVP , (9b) 

6 igee i 

Arrgeea geevP'^grr' 2,2, 
Otgee = ^ 9ee + P gee - i^geevp , (9c) 

grr ogrr O O 

c, . 2aArr'^ vP^'grr'Arr 2vP''gee'Arr 2 , , 2ax{grr'f 

OfArr ~ KaArr r -Vp Arr + 2P Arr + 



3grr 3gee 3 3^^ 

axigee'f a{x'f 2grrax , a ' , j^n "Xffrr'ffee' , xgrr'a' 

5 2 7. 7. V P Arr + ^^rr^Xr h — 

ogee DX ogee o 2grrgee ogrr 

^ Xgee'a' _ agrr'x' _ agee'x' _ 2a'x' _ QXffrr" _^ axgee" _ 2xa" ^ ax" 

'igee 6grr 6gee 3 Sgrr 3gee 3 3 

3aArr^ , K'^a , xgrr'a' xgee'a' a'x' X"" 



rr 



(9d) 



OtK = ^ + + p A + — h , (9e) 

Zgrr O ^grr grrgee ^grr grr 



„ -pr _ ^P'' {gee'f , Arragee' vP'''gee' , P''' gee' , ^^T^r' , ^rrag-rr' ^aK' 2Arra' 

Otl — 5 1 ^ 1 h P i H ^ 

grrgee grr gee ogrrgee grrgee grr^ ogrr grr 

Vgrrfjr grr'P"" SArrax' , vP'-grr" , vP'-geo" , , P""' 

^ — o;r^ 9 — 2 2 — ~R — 2 — Q Q ' ' 

Zgrr ^grr grr X ^grr ogrrgee ogrr grr 



Primes denote derivatives with respect to r. Note that I am following the common practice of using on the 
right-hand side only if it appears differentiated. 

The Hamiltonian constraint is defined by 7i = — KabK"^^ + R and the momentum constraint is defined by 
A4a = DfjK^ — DaK. (Indices on Kab are raised with the physical metric.) With BSSN we also have constraints 
that arise from the definition of the conformal connection functions: Q"^ = T"' — g^'^^bc- -'^^ spherical symmetry these 
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2gr 



constraints become 

Mr = 



2K^ 



5x 



/2 



2x" , 2x 2x9ee" , 2x'50e' , X9rr'9ee' x'9rr' 



X9BS 



/2 



2X5rr 5rr 399 9rr9e0 



A.j^iy '2iI'C 3j4^^^ S-zl-^-p^^^ j^j-^^^^P 



9rr 
9r 



29r 



9ee' 
9rr9ee 



'^9rr9e9 



and the constraint evolution system is 

dtU = p-n' + \aKn - '^^^gr' 

2a 



2ax 



Mr' 



9r 



-Mr . 



9rr9ee 



9rr^9ee 



9r 



29rrgee'^ 



9rr 

I 



ax9rr' 4a'x 2ax9ee' 



9r 



9r 



9rr9ee 



2a'x aX9ee' 
3 3 gee 



Mr 



(10a) 
(10b) 
(10c) 

(11a) 
(lib) 
(He) 



In the three-dimensional application of BSSN one 
must face the possibility that under numerical evolu- 
tion Aab will not remain trace-free. In most three- 
dimensional codes the components Aab are adjusted at 
the end of each timestep to maintain the constraint 
-^ab9°^^ = 0. In spherical symmetry A^b is automatically 
traceless by virtue of the ansatz Eq. Q. As discussed 
previously, the determinant of the conformal metric is 
not constrained in the generalized BSSN formulation. 



IV. HYPERBOLICITY 

Consider the Eulerian (5j^.g = 0) and Lagrangian 
(c?_l5 = — 2D(j/3°) formulations with 1+log slicing Eq. ([l]) 
and either F-driver shift or vanishing shift. With the 
1+log slicing and F-driver shift conditions we can either 
include or exclude the advection terms. These cases will 
be labeled E for Eulerian and L for Lagrangian. A su- 
perscript -|- or — will denote the inclusion or exclusion 
of the advection term in the l-|-log slicing condition. A 
subscript -I- or — will denote the inclusion or exclusion of 
all advection terms in the F-driver shift condition. For 
example, E'\_ represents the Eulerian case with l-|-log slic- 
ing that includes the advection term and F-driver shift 
that includes all advection terms. Similarly, denotes 
the Lagrangian case with 1+log slicing that includes the 
advection term and F-driver shift that excludes all ad- 
vection terms. I also use a subscript to denote vanishing 
shift. In these cases the advection term in the 1+log slic- 
ing condition vanishes, so the superscript can be omitted. 
Thus, Eq and Lq refer to the Eulerian and Lagrangian 
cases with vanishing shift. 

The results on hyperbolicity reported here were ob- 
tained using pseudo-differential techniques [Ml [19]. 
These techniques apply to systems of PDE's that have 
first order time derivatives and second order space deriva- 



tives. Pseudo-differential methods can only distinguish 
between weak and strong hyperbolicity; other methods 
must be used to check for symmetric hyperbolicity. 

For the cases Eq, Lq, which have vanishing shift, the 
full system (equations of motion and 1+log slicing) is 
strongly hyperbolic. For the cases with nonvanishing 
shift, the full system (equations of motion plus gauge 
conditions) is strongly hyperbolic as long as certain in- 
equalities among the field variables are satisfied. 

There are three characteristic fields that are common 
to all cases. These are: 



F'^T 



2\/5rrX 
1 

"^grrgoe 

2 , 

— X ~ 

X9rr 



gee' ± 



25rrX' 

1 



2ff; 



\/9rrX 



zK 



-gee 



(12a) 
(12b) 



The fields ( 12 1) have proper speeds ±1 as measured by 



observers at rest in the spacelike slices. The field ( 12 d) 
has vanishing speed. All cases in which the advection 
term is included in the 1+log slicing condition (or the 
shift vanishes) have characteristic fields 



a' ± ^2agrrlx'i^ 



(13) 



with proper speeds ±-\/2/a. 

The remaining characteristic fields and speeds depend 
on the details of the formulation and the gauge conditons. 
For Eq there are two characteristic fields in addition to 
those displayed in Eqs. (12 1, (13 1. Both are zero speed 
modes. For Lq, the remaining characteristic fields have 
speeds and ff , where = \/ 9rr I xP^ I is the proper 
length shift per unit proper time. 

For the cases E'^ and there are four character- 
istic fields in addition to those displayed above. For 
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the extra characteristic fields have speeds 0, 0, 
±-y/3/(a^x) /2. This system is strongly hyperboHc as 
long as 8ax ^ 3. For L^, the remaining characteristic 
fields have speeds 0, ±y/l/{a'^x) and the require- 
ments for strong hyperbohcity are ^ygrr\fi^\ 7^ 1 and 
2ax ^ 1. For the case E^, the remaining characteris- 
tic fields have speeds 0, /3^ [/3'' ± 0"-)^ + 3/{a^x)] /2 
and the requirement for strong hyperbohcity is {Sax ~ 
3) 7^ ±4-^2ax(7rr/?''- For the case L^, the remain- 
ing characteristic fields have speeds /3'', 15"^ , /?'" ± 

(Z?'")^ + 4/(0:2%) ^2 and the requirement for strong hy- 
perbohcity is 2ax— 1 ^ ±y/2axgrr(3^ ■ In all cases, strong 
hyperbohcity can be spoiled if one or more of the fields 
<?rr, gee, X or a vanishes. I have not carried out a com- 
plete analysis of the restrictions on hyperbohcity for the 
cases Eg and Lg. 

The restrictions on strong hyperbohcity are not a con- 
sequence of the reduction to spherical symmetry. Pre- 
cisely the same restrictions are found for the full three- 
dimensional GBSSN systems. In three dimensions one 
analyzes hyperbohcity by choosing a unit normal cov- 
ector Ua and projecting the principal parts of the equa- 
tions of motion in directions normal and tangential to n^. 
In this way the principal symbol splits into blocks that 
transform as scalars, vectors, and trace-free tensors un- 
der rotations about the normal direction. The vector and 
tensor blocks are strongly hyperbolic without restriction. 
The scalar block in three dimensions is identical to the 
principal symbol for the spherically symmetric GBSSN 
system. It follows that GBSSN has the same restrictions 
on strong hyperbohcity in three dimensions as it has in 
spherical symmetry. 

It should be noted that for traditional BSSN with 
1+log slicing and F-driver shift, including all advection 
terms, strong hyperbohcity requires 2ax 7^ 1 P,9 . For 
traditional BSSN with 1+log slicing and F-driver shift 
and various combinations of advection terms, the restric- 
tions on hyperbohcity have been analyzed in Ref . [20 . 

It is worth noting that the speeds for some of the 
characteristic fields can exceed the speed of light. For 
example, the system E^ has a mode with speed \/2/a 
that is superluminal for a < 2, and a mode with speed 
y/3/{a'^x)/'^ that is superluminal for a^x < 3/4. These 
modes are associated with the 1-l-log and F-driver gauge 
conditions [ST. As such, they appear in any formula- 
tion of the Einstein evolution equations that uses l-|-log 
slicing and F-driver shift. 



The constraint evolution system Eq. (11) is strongly 
hyperbolic with characteristic speeds 0, ±1. The respec- 
tive characteristic fields are 

W + X^^ (14a) 

n T Wx/grrMr + ^xG^^' ■ (14b) 

In three dimensions, the constraint propagation systems 
for both the GBSSN equations and the traditional BSSN 



equations are strongly hyperbolic with causal character- 
istic speeds and ±1. 



V. BOUNDARY CONDITIONS 

The one-dimensional code described in this paper 
could be used to evolve a single black hole with exci- 
sion, or to evolve smooth spherically symmetric fields in 
, provided appropriate boundary conditions are im- 
posed. For the case of fields that are smooth at the ori- 
gin, such as the fields that describe a spherically sym- 
metric star, boundary conditions have been discussed in 
detail in Refs. and elsewhere. In this paper I will 

focus on boundary conditions that allow for black hole 
puncture evolution. In this case the origin r = defines 
the puncture and is also a coordinate singularity in the 
underlying spherical coordinate system. The grid is cell 
centered with spacing Ar. The grid points are located at 
coordinate radii r{j) = {j — 1/2) Ar, where j — 1,2, . . .. 

In applications involving smooth fields with spherical 
coordinates one would normally impose smoothness con- 
ditions at the origin [53]. Specifically, scalars and 
type 2 tensors would be refiection symmetric and vectors 
would be antisymmetric. With puncture evolution we 
need to allow for the possibility that the fields will not 
remain smooth, even if they are smooth initially. We can 
nevertheless attempt to impose smoothness conditions at 
the puncture. In my experience this leads to an unstable 
code. 

To obtain appropriate boundary conditions for punc- 
ture evolution I take a minimalist approach. That is, 
my goal is to impose as little information as possible 
through boundary conditions so that the puncture can 
evolve without unnecessary infiuence from the outside. 
After some experimentation, I found that the following 
prescription works well for the five Lagrangian cases Lq 
and L^. For the variables grrj Xj 

K, F'', a, and B'', 

no boundary conditions are imposed. That is, for grid 
points near the puncture, the finite difference stencil is 
shifted away from the puncture so that no guard cells are 
needed. 

Typically, the condition I impose on ggg and /?'' is that 
these variables should vanish at the puncture. In keep- 
ing with the minimalist approach I only use one layer of 
guard cells. The guard cells {j = 0) are filled via the 
relations 



geeiO) 



h315<?ee(l) + 210gggi2) - I26gegi3) 
+455ee(4)-75ee(5)]/63 , 
[-315/3''(1) + 210/3''(2) - 126/3''(3) 
+45/3'^ (4) -7/3'- (5)] /63 , 



(15a) 



(15b) 



where the numbers in parentheses denote grid points. 
These relations insure that ggg and f3^ vanish at the punc- 
ture to sixth order accuracy. In other words, a fifth order 
polynomial fit across the grid points j — . . .4 with value 
( 15 1 at j = will yield a function that vanishes at r = 0. 
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FIG. 1: The proper speed \/ {P^Y'Qrr Ixl of the coordinate 
system, plotted as a function of coordinate radius r at times 
t = M, 2M, 3M, and 4Af. 



The results presented in Sec. VII are obtained from a 
code that uses fourth order accurate spatial differencing. 
In the bulk of the computational domain, I use centered 
stencils that extend across five grid points for first and 
second spatial derivatives. (An exception to this rule is 
made for advection terms like /3'"9r5rr- For these terms, 
the finite difference stencil is shifted by one grid point 
in the "upwind" direction.) Since only one guard cell is 
available at the boundary r = 0, the stencils near r = 
are shifted so that they require just one grid point on the 
side closest to the puncture. 

For the Eulerian cases Eq, and E^, no boundary 
conditions are imposed on the variables gm gee, X, ^rr, 
K, r*", a, and _B''. The only variable that remains is the 
shift vector /?''. For the shift I impose drP^ — via the 
relationship 

/3-(0) = [17/3'-(l) + 9/3'-(2) - 5/3^3) + /3''(4)]/22 . (16) 



Here again the numbers in parentheses label grid points. 
A fourth order polynomial fit across the grid points j = 
. . . 4, with value ( 16 ) at j = 0, will yield a function with 
vanishing derivative at r — 0. 

For the cases E^ and EZ with the condition Eq. (16 1, 



an instability develops at the puncture within a time of 
a few M. I have not found suitable inner boundary con- 
ditions for these cases. 

The minimalist approach to boundary conditions at 
r = is not unreasonable. Simulations show that as the 
data evolve from the initial conditions, the light cones at 
grid points near the puncture quickly tip outward (to- 
ward r — 0) with respect to the time flow vector field. In 
Fig. ([ij the proper speed of the coordinate system with 
respect to the spacelike slices is plotted at various times 
up to t = 4M. During the first 4M of evolution the hori- 
zon moves from r = M/2 to just beyond r = M. From 
the graph we see that hy t = M all of the grid points 
inside r ^ M/3 are moving faster than the speed of light 



(which is unity). By t = AM the region of superlumi- 
nal grid movement has expanded to about r k M, the 
location of the horizon. There is little change beyond 
t = AM. 

Once the time flow direction moves outside the light 
cone, one might expect that any boundary conditions 
imposed at the puncture would be irrelevant. This is 
not quite correct because GBSSN (as well as traditional 
BSSN) with l-l-log slicing and F-driver shift contains 
characteristic modes that can travel faster than light 
For example, the 

±y27a and ±^3/(4a2x) 



E^ system has modes with speeds 
The system has modes 

with speeds -^^J^Ja, ±-\/l/(a^x)i ^'^^ ■ 

Ideally, one should impose boundary conditions on all 
incoming characteristic fields at the puncture. Whether 
or not a characteristic field is incoming or outgoing de- 
pends on the proper speed of the mode relative to the 
proper speed of the coordinate grid, which is 13^ . For a 
typical puncture evolution the modes listed in Eq. (13) 
with speeds -^^J^Ja are outgoing, at least for times be- 
yond the first few M. For £^t, the characteristic field 



^figr 



8ax 



4a 



^/3g^{8ax — 3) Sax ~ 3 



K (17) 



has speed ±\/3/(4a^x) ^^^d is incoming at the punc- 
ture throughout the evolution. Thus, we would expect 
that certain boundary information must be imposed at 
the puncture to fix this incoming mode. The condition 
Eq. ( 16 1 is not ideal, but appears to be sufficient to allow 
for stable evolutions. The case is more difficult to 
analyze, because in addition to an incoming mode it also 
contains a mode that travels along the puncture bound- 
ary. 

The above discussion shows that for puncture evolu- 
tion in one dimension, some boundary conditions at the 
puncture are necessary. This might come as a surprise 
since, apparently, no boundary conditions are imposed at 
the puncture in three-dimensional codes. This is proba- 
bly an illusion. More likely, the incoming modes in the 
three-dimensional case are fixed implicitly through the 
numerical (finite differencing) scheme used in the vicin- 
ity of the puncture. 

The details of the finite differencing scheme in the 
vicinity of the puncture would have little affect on the 
data away from the puncture, and would be especially 
difficult to detect outside the black hole. Evidence for 
this assertion can be found in recent work on the "tur- 
ducken" approach to black hole evolution [T^ [23], 
in which the black hole interior is stuffed with artificial, 
constraint- violating data. In this case it is found, both 
for one-dimensional simulations with E^ and for three- 
dimensional simulations with traditional BSSN, that the 
stuffing can only affect the slicing and coordinate con- 
ditions outside the black hole [15]. The slicing can be 
changed outside the black hole d ue t o the presence of the 
mode Eq. (13) whose speed •\/2/a can become super- 
luminal. However, inside a coordinate radius of about 
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r ~ 0.2M, this mode moves more slowly than the coordi- 
nate grid and cannot propagate to the black hole exterior. 
The coordinate conditions outside the black hole can be 
affected by the interior data through the presence of the 
mode Eq. (17). The conclusion is the following: the de- 



tails of how the puncture is handled should not affect the 
spatial geometry or slicing of the black hole. Therefore, 
the boundary conditions at the puncture, whether they 
are imposed explicitly in one dimension or implicitly in 
three dimensions, do not affect the physics. Their only 
effect is to change the details of how the spatial coordi- 
nates are shifted in the spacelike slices. 

Another subtle issue concerning the origin r = is how 
the variable F'' is treated. Because F'" is singular at the 
puncture, F*" ~ — 2/r, one would expect the finite dif- 
ference calculation of F*"' to generate large errors. This 
is indeed the case, but the errors are effectively trapped 
near the puncture. 1 have carried out numerical simula- 
tions using the "regularized" variable F^^^ = F'' -I- 2/r in 
place of F''; a similar approach is advocated in Ref. 
In some cases this does help reduce the errors near the 
puncture, but it does not change the accuracy of the code 
elsewhere. 

Finally, let me mention that at the outer boundary of 
the computational domain I use two layers of guard cells 
and simply freeze the data at those points. In general the 
outer boundary is placed at a coordinate distance that is 
greater than the total run time. Assuming causal prop- 
agation, the outer boundary conditions do not affect the 
fields near the puncture during the course of a simulation. 



VI. KRUSKAL-SZEKERES DIAGRAM 



Note that f3r = P^grr/X- Let 14 = const be a null curve 
in the r-t plane. By solving dL{ = dtU dt -\- drU dr — 
for dt and inserting the result into ds^ = 0, we can 
easily show that U satisfies d±U — ztdrU. Here, the 
operator d± = {dt — (3'^dr)/a is the normalized deriva- 
tive orthogonal to the spacelike hypersurface acting on 
scalars. The operator dr = y/ xl 9rrdr is the normalized 
derivative tangent to the spacelike hypersurface. Now 
observe that u ± v are null coordinates. It follows that 
u and V satisfy the equations d±{u + v) = dr{u + v) and 
— v) — —dr{u — v), where the signs are chosen so 
that drU is positive. 

The reasoning above shows that the Kruskal-Szekeres 
coordinates satisfy the advection equations 



dtiu + v) 
dt{u - v) 



(/?'■ + a^xl9rT)dr{,u -t- v) , (21a) 
iP'- -a^x/9rr)dr{u-v) . (21b) 



We can integrate these equations forward in time along 
with the BSSN variables. To do so we must choose ini- 
tial data and boundary conditions. The initial data for 
puncture evolution of a single black hole is obtained from 
the Schwarzschild metric in isotropic coordinates. 



9rr 

gee 



= 1 , 



X = (l + M/2r)- 



(22a) 
(22b) 
(22c) 



along with vanishing extrinsic curvature. The areal ra- 
dius is i? = \fgeeTx = ?'(1 + M/2rY. If we assume that 
the initial slice is w = 0, Eq. (19 1 shows that u = F{r) 
where F{r) is defined by 



One of the benefits of assuming spherical symmetry 
is that we can visualize the motion of the slices and 
grid points on a Kruskal-Szekeres diagram. In Kruskal- 
Szekeres coordinates, the Schwarzschild black hole metric 
is [251 



ds^ 



32M3 
R 



{~dv^ + du^) + dn^ . (18) 



Each point in the u-v plane is a sphere of areal radius R, 
where R is defined by 



\2M 



(19) 



Our goal is to track the numerical grid points as they 
move through the u-v plane. 

There are a number of ways one can relate the numer- 
ical data to the u-v values of the grid points. 1 have had 
the most success with the following approach. Consider 
the metric for a spherically symmetric spacetime in 3+1 
notation: 



ds^ = -{a^ - (3'-f3r)dt^ + 2l3rdtdr + ^dr^ 

X 



gee 



dn^ 



(20) 



F{r) 



2M 



1 



M 
2^ 



„r(l+A//2r)V(4M) 



(23) 



More generally, we can take the initial slice to be boosted 
in the Kruskal-Szekeres diagram: 



u = F{r) cosh {to /AM) 
V = F(r) sinh (<o/4Af ) 



(24a) 
(24b) 



The constant tg represents the invariance of the 
Schwarzschild geometry under the action of the Killing 
vector field vd^ + udy. If we take ^0, then the initial 
data slice is related to ?; = by sliding the slice along the 
orbits of the Killing vector field. 

The freedom to choose Iq is useful, because starting 
from u = the slice quickly becomes stretched around 
the singularity i? = 0. That is, on a Kruskal-Szekeres di- 
agram, the slice becomes visually difficult to distinguish 
from the singularity after a few M of evolution time. 
We can avoid this somewhat by starting with an initial 
slice with < 0. In that case the early time slices are 
difficult to visualize, but the slices at times i « tg are 
clearly displayed in relation to the singularity and hori- 
zon. Another useful technique is to slide the slice along 
the Killing vector field after each timestep. This can be 
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used, for example, to keep grid points from crossing the 
V axis, or to fix the location in the Kruskal-Szekeres di- 
agram where the slices cross the horizon. 

The variables u±v are left or right -moving, depending 
on the sign of f3^ ± a^/x/g^. For puncture evolutions, 
within a coordinate time of a few M, the shift vector 
grows and the lapse collapses so that dominates over 
a\/xl 9rr near r = 0. Therefore, beyond a time of a 
few M, both characteristic fields are outgoing (toward 
decreasing r) at the boundary r = 0. In the vicinity 
of the origin I use one-sided differencing for u and f , 
appropriate for outgoing fields. Thus, no guard cells are 
needed for u and w at r = 0. 

At the outer boundary of the computational domain 
we have w and xl 9rr ~ 1- Thus, the char- 
acteristic field u + V IS incoming and the field u — v is 
outgoing. Boundary conditions are applied with the fol- 
lowing scheme. First, let us equate the coefficients of dt^ 
in Eqs. (18) and (20 1 to obtain 



We also find, by differentiating the definition (19) 



R 



(25) 



32M 

where the dot denotes dt- Now differentiate the definition 



( 19 1 with respect to time to obtain a second equation 
involving ii and v. Together these equations have the 
solution 



u = ^[uG + ev^G-^ + ^GHia^-p-^Pr)] , (26a) 



V = —\^G + eu^G-^ + ^GH{a-^~l3-l3r)J ,(26b) 

where G = {R/2M - 1) exp(i?/2M), H = 
(i?/32M3)exp(i?/2M), and e = u/\u\. Note that 
G can be written in terms of R and R. Also recall that 
the areal radius can be written in terms of the BSSN 
variables as i? = \/gee/X) so the time derivative R 
depends on dtgge and dtX- The time derivatives of these 
BSSN variables can be replaced with the right-hand 



sides of the BSSN equations ( 9a I and ( 9c 



Equations ( 26 1 are a set of ordinary differential equa- 
tions for u and v. These can be integrated forward in 
time without the calculation of any spatial derivatives. I 
use these equations to determine u and v in the vicinity 
of the outer boundary, providing a layer of two guard 



cells for the interior calculation based on Eqs. (21a 



One can try to use Eqs. (26) to evolve u and v every- 



where. In practice this does not work well. It is difficult 
to integrate these equations across the horizons where 
G = 0. Also it is difficult to maintain accuracy and sta- 
bility as grid points evolve across the v axis, where the 
value of e changes from —1 to +1. 

Finally, let me mention a method for placing the data 
of a single time slice into a Kruskal-Szekeres diagram. 
From the numerical data one can compute the areal ra- 
dius R and the proper distance from the horizon, L, for 
each grid point. In the u-v plane the metric gives 



32AP 
R 



-R./2M 



(Au2 _ At;^) 



(27) 



uAi 



vAv = 



R 
8M2 



(28) 



These equations can be solved for Au and Av as functions 
of AR and AL. The values of u and v for each grid 
point are found by numerical integration, starting from 
some initial values for u and v. These initial values must 



satisfy the restriction (19), which leaves some freedom of 



choice. For example, if the intitial areal radius is less than 
2M, one can choose the initial values to be u = and 
V = -^/l — R/2Me^/'^^^ . The freedom in choosing the 
initial values for u and t; is a consequence of the Killing 
isometry of Schwarzschild spacetime. 



VII. RESULTS 

There are a large number of cases that we can consider, 
depending on our choice of Eulerian versus Lagrangian 
evolution for the advection terms that we include or 
omit in the slicing and coordinate conditions, and the 
value for the damping parameter 77 that appears in the F- 
driver shift. I will not attempt to provide a complete ex- 
amination and comparison of all these possibilities here. 
Rather, I will give a sampling of some of the interesting 
results. 

The initial metric for puncture evolution is given in 
Eq. ( 22 1 . Initially, the conformal connection function is 
F'' = — 2/r, and the extrinsic curvature components A„ 



and K vanish. The shift vector fi"^ and the auxiliary field 
are chosen to vanish at the initial time. The initial 
lapse function is either unity, a = 1, or has the "pre- 
coUapsed" form a — {\ + M/2r)~2. In all cases the code 
uses fourth order finite differencing in space and fourth 
order Runge-Kutta for time stepping. 

The geometrical description of puncture evolution for 
single black holes is now well understood [T31 [HI US) [271 
,28] . The initial geometry is a wormhole. As the evolution 
begins, the grid points near the "other" asymptotically 
flat end quickly shift into the black hole interior. The 
numerical data then settle onto a portion of a "trumpet 
slice" of the black hole. Such a slice asymptotes to a 
constant areal radius of about 1.3M. The key ingredient 
responsible for this behavior is the F-driver shift condi- 
tion. Tests using the one-dimensional code with F-driver 
shift invariably show the same qualitative evolution for 
the physical geometry, the same evolution from wormhole 
to trumpet, for all of the cases and L^. However, the 
detailed behavior of the GBSSN variables can differ sig- 
nificantly from case to case. 

Figures ([2||5| show the metric components grr, gee, 
the conformal factor Xi E^nd the shift vector at time 
t = 50M for the cases i?^ and L\. The data for these 
graphs were taken from simulations with resolution Ar — 
Af/100. The differences near the puncture are evident. 

One might guess that the difference between the cases 
Et and L\ is a consequence of the different boundary 
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FIG. 2: The conformal metric component Qrr near the punc- 
ture r = at time t — 50M for the systems Ef and Lt. 




FIG. 5: The shift vector P'' a.t t = 50M for E+ and L+. 




FIG. 3: The conformal metric comonent 
E+ and L+. 



at t = 50M for 



X 




FIG. 4; The conformal factor x at t = 50M for E^ and L:]; 



conditions imposed at r = 0. However, the Eulerian 



equations cannot be evolved with the Lagrangian bound- 
ary conditions (15 1. Recall that the Lagrangian bound- 
ary conditions are designed to keep gee and /?'' equal to 
zero at the puncture. If these conditions are used with 
the Eulerian evolution equations, the values of gee and 
f3^ near the puncture grow rapidly and a sharp gradient 
develops between the first few grid points j = 1,2, .. . 
and the guard cell j — 0- The code crashes shortly after 
t = lAI. With the boundary conditions (16 1, the fields 



gee and /?'" are allowed to develop nonzero values at the 
puncture. Figures (|3| and (|5| show that indeed gee and 
Z?*" evolve to nonzero values in the Eulerian case E^. 

The converse is also true: The Lagrangian equations 
will not evolve stably with the Eulerian boundary condi- 
tions (16 1. If one attempts such an evolution, the code 



quickly develops problems at the puncture and crashes. 
It appears that the behavior of the fields near the punc- 
ture is primarily dictated by the equations of motion in 
the bulk. Perhaps the main role played by the boundary 
conditions ( 15 1 or ( 16 1 is to help insure numerical stabil- 



ity for the GBSSN equations in spherical symmetry. 

Although the results obtained with and ap- 
pear quite different, they represent the same slicing of 
the same spacetime geometry. It is clear that the space- 
time geometry is that of a Schwarzschild black hole, since 
a Schwarzschild black hole is the unique solution of the 
vacuum Einstein equations with the chosen initial data. 
It is not quite so obvious that the evolutions E^ and 
lead to the same foliation of that spacetime geome- 
try. The difference between the Eulerian and Lagrangian 
cases resides in the choice for dj_ \ng, as seen in Eqs. (|5|. 
The terms containing d± In g can be absorbed into redef- 
initions of the field variables. For example, Eq. ([sji) can 
be written as d^cf) — —aK/6 where (f> = (f>+ln g^^^'^ . Sim- 
ilarly, Eq. pp) can be written as d±gnh = ~2aAab where 



<jab 



9 



gab- The physical metric ([3p) is invariant un- 

This invariance is 
a consequence of the conformal symmetry of the GBSSN 
equations (16.j. In this way we see that, independent of 



der this redefinition: e^'^gab — e'^'^gab 
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FIG. 6: Penrose diagram showing the grid points for the cases 
(large green circles) and (small red circles) at time 
t = 2.5M. The data for blend together to form a thick 
curve that ends in the middle of the diagram. 




10 20 30 40 50 60 
r/M 

FIG. 7: grr as a function of coordinate radius r at time t = 
50M for the case L1. The three curves correspond to values 
0, 1/M, and 2/M for the damping parameter rj. The peak 
values increase with increasing 77. 



the choice for d± In g, the physical metric and extrinsic 
curvature undergo identical evolutionary paths. 

Figure ^ shows a penrose diagram with the grid 
points at time t — 2.5M as obtained from the and 
cases. The slice is the same for the two cases. What is 
clearly different is the coverage of that slice by the avail- 
able grid points. In the Eulerian case, the relatively large 
positive shift near the puncture, seen in Fig. ([5|, drives 
the grid points away from the left spacelike infinity and 
into the interior of the black hole. In the Lagrangian case 
the grid points are also driven from the left spacelike in- 
finity into the black hole interior, but not as quickly. 

The behavior of the three-dimensional GBSSN equa- 
tions in the case appears to be close to the behavior 
of three-dimensional BSSN codes. It is particularly use- 
ful to compare the results presented below with those of 
Ref. [7 . In particular. Figs. ([7 10 1 show graphs of the 
variables grr, Xi ™d as functions of the coordinate 
radius r at time t = 50M. The evolution is type (La- 
grangian case with 1-l-log slicing including the advection 
term and F-driver shift excluding advection terms), and 
begun with a pre-coUapsed lapse. The resolution used 
for these simulations was Ar = M/100. The curves 
in these figures show the results obtained with three val- 
ues of the damping parameter, rj = 0, 77 = 1/M, and 
T] = 2/M. 

Recall that rj appears in the F-driver shift condition. 
For l-|-log slicing including the advection term the shift 
vector does not affect the slicing of spacetime or the spa- 
tial geometry, it only changes the coordinate system (or 
distribution of grid points) from one slice to the next 
|13j . Thus the results obtained with different values of 
r] are physically equivalent. Nevertheless, there is rea- 
son to prefer a value that will allow the lapse and shift 
to adjust to the Killing symmetry of the Schwarzschild 
geometry. That is, we want the gauge to be "symmetry 




FIG. 8: X a-s a function of r at time t = 50Af for the case L_ 
with rj — 0, 1/M, and 2/M. The peak values increase with 
increasing rj. 



seeking" and bring the data close to a stationary state at 
late times [21]. For the one-dimensional code with 1-1- log 
slicing the value 77 = is a good choice. With ry = the 
data settle very close to a stationary state by t = 50M. 

With positive 77 the data continue to evolve indefinitely. 
This tendency has been seen in three-dimensional sim- 
ulations [7| where it is described as "coordinate drift". 
We can monitor this drift by computing the L2 norm of 
the right-hand sides of the GBSSN equations. (The L2 
norm is defined as a sum over the right-hand sides of the 
evolution equations for grr, X, , ct, ^'^d f3^, and a 

sum over grid points between r — and r = 25M.) For 
each value of 77 the L2 norm rises to a peak value just be- 
yond t = 2M. For 77 = the norm drops slowly between 
t « 2M and t w 25, then drops rapidly for t > 25M. 
The sudden change in slope at i « 25M occurs for the 
following reason. At early times there is a relatively large 
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FIG. 9: Q as a function of r at time t = 50M for the case L_ 
with Tj — 0, 1/M, and 2/M. The peak values decrease with 
increasing 77. 



FIG. 11: Common Logarithm of the L2 norm of the right- 
hand sides of the equations of motion versus time. 
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FIG. 10: as a function of r at time t = 50M for the case 
L1 with ?7 = 0, 1/M, and 2/M. The peak values increase 
with increasing rj. 
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FIG. 12: Convergence plot for x at t = 2Af. The data are 
obtained from simulations at resolutions Ar = Af/50, M/lOO, 
M/200, and A//400. The curves are scaled by powers of 16, 
appropriate for fourth-order convergence. 



adjustment in the geometry near the puncture as the grid 
points near the puncture pull away from spacelike infin- 
ity, enter the black hole interior, and relax to the "trum- 
pet slice" [13]. This adjustment can be seen in the values 
of the right-hand sides as a pulse that propagates out- 
ward. At t ~ 25M the pulse passes beyond the region 
in which the L2 norm is computed. The rapid decay for 
77 = beyond t « 25M shows that the data in the in- 
terior region r < 25M of the computational grid is very 
effectively becoming stationary. 

The code described here uses fourth order finite differ- 
encing in space and time. However, the code does not 
always exhibit fourth order convergence near the punc- 
ture. Figures (12 1 and (13 1 show convergence plots for 
X at early and late times, t — 2M and t = 50M, for the 
case L^. Each graph shows three curves, obtained by 
computing differences between values of x at successive 
resolutions and scaling by powers of 16. To be explicit, let 



us add a subscript to x to denote the resolution Ar. Then 
the curves in Figs. ^ and ^ are (xm/50 - Xm/wo), 
16(xm/ioo - XM/20o)7and 256(xm/200 - Xm/40o)- These 
graphs show that at early times the data are fourth-order 
convergent, while at late times the data are fourth order 
convergent outside the region close to the puncture. 



Figure ( 14 1 is a convergence plot for x at the interme- 
diate time t = 15M, again for the system L^. This graph 
shows four curves, {xm/50-Xm/ioo), 4(xm/ioo-Xm/20o), 
16(xa//200 - XM/400), and 64(xm/4oo - Xm/soo)- The 
curves are scaled by powers of 4, which is the scaling 
expected for second-order convergence. An error pulse 
is generated near the puncture and propagates outward 
through the computational domain. Note that this error 
is delayed at higher resolution, resulting in a shift in the 
curves. 

These tests suggest that the puncture generates an ap- 
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FIG. 13: Convergence plots for x&tt = 50Af. The resolutions 
and scaling are the same as those used in Fig. ( 12 1. 
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FIG. 14: Convergence plot for x at t = 15M. The resolutions 
are the same as those used in the previous figures. In this 
case the curves are scaled by powers of 4, appropriate for 
second-order convergence. 



proximately second-order error of the form 

XAr = Xcxact + £{r - t + 2A-/ log(A//Ar)) (29) 

r 

for some function S. The pulse appears to form from 
errors near the puncture. Note that the horizon begins 
at coordinate radius r = A//2 and stays within r < M 
throughout the evolution. The pulse moves with a speed 
of approximately 1. It has an amplitude approximately 
proportional to Ar^/r. At higher resolution the pulse is 
delayed. This delay appears as a shift toward the punc- 



ture in Fig. ( 14 1. The amount of delay or shift is propor- 



tional to the logarithm of 1/Ar. 

Similar convergence tests with the system do not 
show this anomolous behavior. This suggests that the 
origin of the non-fourth order error for L'^ is related to 

the presence of modes with the unusual proper speed /J"". 



Because the characteristic fields and speed for traditional 
BSSN are more closely related to those for E^, I would 
not expect to see such behavior in current 3D codes using 
the traditional BSSN equations. 



Figure ( 15 1 shows a plot of the L2 norm of the con- 
straints (10 1 as functions of time. The constraint val- 
ues are stable over time, with only the Hamiltonian 
constraint Ti. showing a slight upward drift. All of the 
constraints exhibit large non-convergent errors near the 
puncture. These errors dominate the calculation of the 
L2 norm over the entire computational domain. As a 
result, the data for Fig. (151 was obtained by omitting 
the region < r < 0.5M. The data beyond r > 25M 
is also excluded so that outer boundary effects are ig- 
nored. The data for the momentum constraint and 
the conformal connection function constraint show 
nearly perfect fourth order convergence over this region 
M/2 < r < 25M. The results for the Hamiltonian con- 
straint H are not so ideal. I suspect this is related to 
errors that arise in computing Ti., which depends on sec- 
ond spatial derivatives of the metric components ggg and 
X- The constraints Mr and do not involve second 
order derivatives. 

The flow of grid points on a Kruskal-Szekeres dia- 
gram can be particularly enlightening. Figure ( 16 1 shows 



the grid points at t — 50M located on a portion of the 
Kruskal-Szekeres diagram. This data was obtained from 
a simulation with rj = 0, and initial lapse a = 1. The 
initial slice in the diagram was obtained from Eq. (24 1 
with ^0 = — 50Af . The resolution for this simulation was 
Ar — Af/100. For clarity, only the odd numbered grid 
points are shown in this figure. The smooth curve was 
obtained from a simulation with vanishing shift vector, 
case Lq. 



The first grid point in Fig. (16 1, j = 1, is at location 
« —3.7, V ~ 3.7 between the horizon and the physical 
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FIG. 15: Common log of the L2 norm of the constraints TC, 
Mr, and . For each constraint two resolutions are plotted, 
Ar = Af/100 and Ar = A'//200. For the higher resolution, the 
constraints are multiplied by a factor of 16. The two curves 
for Mr overlap one another, and are nearly indistinguishable 
in this figure. Likewise the two curves for overlap one 
another. 
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singularity. As shown in Ref. |13j . the grid points near 
the puncture are rapidly drawn into the black hole in- 
terior by the F-driver shift condition. This is, however, 
a resolution -dependent effect. In the limit of infinitely 
high resolution, the grid points merge into the smooth 
curve. That curve crosses the horizon and asymptotes to 
spacelike infinity (u ~ — oo, v ~ finite). 

Figure (17) compares the results of simulations with 
cases LZ and Lq. Again, we use 77 = 0, unit initial lapse, 
and choose to = — 50M. The figure shows grid points 
j — 3, 5, 7 . . . at time t = 50M for the LZ evolution. 
The grid point j — 1, which is not shown, is at location 
u « —12.7, V ~ 12.7. The smooth curve on the Kruskal- 
Szekeres diagram was obtained from the Lq evolution. 
Note that in this case the grid points do not lie on the 
curve. This is because without the advection term in the 
1+log slicing condition, the slicing depends on the shift 
vector. 

Finally, let me comment on the results of long term 
evolution. Typically, after a run time on the order of a 
few light crossing times, the code will crash due to the 
development of a shock. For example, with evolution 
on a grid with outer boundary at 125M, steep gradi- 
ents and spikes develop in the BSSN variables at about 
r « 8.5M. The behavior causes the code to crash at 
t « 250M. With the outer boundary placed at 75M, 
the shock develops at r « 6M and the code crashes at 
t w 160M. It is possible that improvements to the outer 
boundary conditions might postpone or even eliminate 



these shocks. On the other hand, it has been suggested 
that gauge shocks are a generic result of l-l-log type slic- 
ing conditions pi . Currently, such behavior has only 
been seen with 1+1 codes. 



Acknowledgments 

I would like to thank Dae-Il Choi, Pablo Laguna, 
Olivier Sarbach and Manuel Tiglio for helpful discussions. 
This work was supported by NSF grant PHY-0600402. 




FIG. 16: Kruskal diagram ai t — 50M for simulations with 
1+log slicing. The hyperbolic curve is the future singularity 
and the lines at +45° are the horizons. 
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